Analyses des VCF en sortie du pipeline
Elodie PARIS
2024-03-17
Import des packages
rm(list=ls())
library('dplyr')
library('tidyr')
library('VennDiagram')
library('grid')
library('futile.logger')
library("ggplot2")
library("gplots")
library('plotly')
library('tools')
library('cowplot')
library('UpSetR')Import des VCF
Import des 3 vcf obtenus après lancement du pipeline (vcf_12x, vcf_24x, vcf40x) et import du vcf de référence (vcf_ref).
vcf_12x = read.table("/home/elodie/Documents/Analysis/gold_12x_on_data_elodie_filtervarianttranches_2D_decomposed_normalized_uniq.vcf")
vcf_24x = read.table("/home/elodie/Documents/Analysis/gold_24x_on_data_elodie_filtervarianttranches_2D_decomposed_normalized_uniq.vcf")
vcf_40x = read.table("/home/elodie/Documents/Analysis/gold_40x_on_data_elodie_filtervarianttranches_2D_decomposed_normalized_uniq.vcf")
vcf_ref = read.table("/home/elodie/Documents/Elodie/Datas_ismael_ref/gold_on_data_elodie_ismael_decomposed_normalize_uniq.vcf")Exploration des VCF et des variants
Nombre et type de variants par VCF
plot1 = ggplot(nombre_variants, aes(fill=type, y=vcf, x=nombre)) +
geom_bar(position = "stack", stat = "identity") +
scale_fill_manual(values = c("#bebada", "#ffffb3", "#8dd3c7", "#170f0f")) +
labs(fill = "Type de variants", x = "Nombre de variants", y = "VCF") +
theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5)) +
labs(title = "Nombre et type de variant par VCF")
ggplotly(plot1)La distribution des types de variants semble équivalente pour chaque VCF. Cependant dans les 3 VCF générés par le pipeline on trouve des variants “Undetermined” mais pas dans le VCF de référence.
Que sont les Undetermined non retrouvés dans le VCF de référence ?
undetermined = vcf_all %>% filter(type == "undetermined")
print(undetermined)## variants chromosome position
## 1 chr4_66388036_A_* chr4 66388036
## 2 chr4_66440407_CATATATATACATATATATACAT_* chr4 66440407
## 3 chr4_66214417_A_* chr4 66214417
## 4 chr4_66440407_CATATATATACATATATATACAT_* chr4 66440407
## 5 chr4_66473213_A_* chr4 66473213
## ref alt format Id gene
## 1 A * 1/.:0,5:6:41:205,41,45 423 EPHA5
## 2 CATATATATACATATATATACAT * 1/.:0,16:22:99:867,243,285 535 EPHA5
## 3 A * 1/.:0,11:15:99:717,152,110 110 EPHA5
## 4 CATATATATACATATATATACAT * 1/.:0,22:31:99:1233,367,435 547 EPHA5
## 5 A * 1/.:0,4:15:99:1163,469,449 651 EPHA5
## type vcf
## 1 undetermined vcf_12x
## 2 undetermined vcf_24x
## 3 undetermined vcf_40x
## 4 undetermined vcf_40x
## 5 undetermined vcf_40xVariants avec * en alt lié au vt décompose ?
pos66388036 = vcf_all %>% filter(variants == "chr4_66388034_GTA_G")
print(pos66388036)## variants chromosome position ref alt format
## 1 chr4_66388034_GTA_G chr4 66388034 GTA G 0/1:2,5:7:46:162,0,46
## 2 chr4_66388034_GTA_G chr4 66388034 GTA G 0/1:7,8:15:99:250,0,267
## 3 chr4_66388034_GTA_G chr4 66388034 GTA G 0/1:13,15:28:99:472,0,494
## 4 chr4_66388034_GTA_G chr4 66388034 GTA G 0/1:.:362:146,111:176,136:155
## Id gene type vcf
## 1 422 EPHA5 deletion vcf_12x
## 2 441 EPHA5 deletion vcf_24x
## 3 453 EPHA5 deletion vcf_40x
## 4 413 EPHA5 deletion vcf_refpos66388036bis = vcf_all %>% filter(variants == "chr4_66388036_A_G")
print(pos66388036bis)## variants chromosome position ref alt format Id
## 1 chr4_66388036_A_G chr4 66388036 A G ./1:0,1:6:41:205,176,204 424
## gene type vcf
## 1 EPHA5 SNV vcf_12xA cette position délétion de 2 bases (66388035+66388036) bien détectée dans les 4 VCF Substitution de A>G en 66388036 dans le VCF_12x (comme peu de reads, variant caller se trompe? Faux positif?)
pos66440407 = vcf_all %>% filter(position == 66440407)
print(pos66440407)## variants chromosome position
## 1 chr4_66440407_CATATATATACATATATATACAT_* chr4 66440407
## 2 chr4_66440407_CATATATATACATATATATACAT_C chr4 66440407
## 3 chr4_66440407_CATATATATACATATATATACAT_* chr4 66440407
## 4 chr4_66440407_CATATATATACATATATATACAT_C chr4 66440407
## ref alt format Id gene
## 1 CATATATATACATATATATACAT * 1/.:0,16:22:99:867,243,285 535 EPHA5
## 2 CATATATATACATATATATACAT C ./1:0,6:22:99:867,631,654 536 EPHA5
## 3 CATATATATACATATATATACAT * 1/.:0,22:31:99:1233,367,435 547 EPHA5
## 4 CATATATATACATATATATACAT C ./1:0,9:31:99:1233,875,889 548 EPHA5
## type vcf
## 1 undetermined vcf_24x
## 2 deletion vcf_24x
## 3 undetermined vcf_40x
## 4 deletion vcf_40xZone avec des reads ayant des dététions différentes qui se chevauchent dans les 4 BAM. Le VCF de référence n’appelle pas de variants dans cette zone. Zone difficile pour le variant caller, faire appel à un variant caller spécialiste des délétions?
pos66214417 = vcf_all %>% filter (position == 66214417)
print(pos66214417)## variants chromosome position ref alt format Id
## 1 chr4_66214417_A_* chr4 66214417 A * 1/.:0,11:15:99:717,152,110 110
## 2 chr4_66214417_A_T chr4 66214417 A T ./1:0,4:15:99:717,468,450 111
## gene type vcf
## 1 EPHA5 undetermined vcf_40x
## 2 EPHA5 SNV vcf_40xpos66214413 = vcf_all %>% filter (position == 66214413)
print(pos66214413)## variants chromosome position ref alt format
## 1 chr4_66214413_TATAA_T chr4 66214413 TATAA T 0/1:1,9:10:15:375,0,15
## 2 chr4_66214413_TATAA_T chr4 66214413 TATAA T 0/1:4,11:15:99:450,0,110
## Id gene type vcf
## 1 108 EPHA5 deletion vcf_24x
## 2 109 EPHA5 deletion vcf_40xDélétion de 4 bases de 66214413 à 66214417 vu dans les 4 BAM mais uniquement appelée dans VCF 24x et 40x. Non appelé dans le VCF de référence. Subtitution en 66214417 aussi un A>T vu dans les 4 BAM mais uniquement appelée dans le VCF 40x. Non appelé dans le VCF référence.
pos66473213 = vcf_all %>% filter (position == 66473213)
print(pos66473213)## variants chromosome position ref alt format Id
## 1 chr4_66473213_A_C chr4 66473213 A C 1/1:0,18:18:54:720,54,0 635
## 2 chr4_66473213_A_* chr4 66473213 A * 1/.:0,4:15:99:1163,469,449 651
## 3 chr4_66473213_A_C chr4 66473213 A C ./1:0,11:15:99:1163,206,135 652
## gene type vcf
## 1 EPHA5 SNV vcf_24x
## 2 EPHA5 undetermined vcf_40x
## 3 EPHA5 SNV vcf_40xDans le VCF de référence en position 66473209 à 66473213 : délétion de 4 bases appelée mais peu de reads avec la délétion et beaucoup plus de reads avec substitutions de A>C en 66473213. Substitutions appelée uniquement dans VCF 24x et 40x. Dans les zones où il y a des délétions et substitutions, le variant caller rencontre des difficultés?
Nombre de variants par gène et par VCF
plot2 = ggplot(nombre_variants_par_gene, aes(fill=gene, y=vcf, x=nombre)) +
geom_bar(position="stack", stat="identity") +
scale_fill_manual(values = c("#fdc086", "#91bfdb", "#b2df8a")) +
labs(fill = "Nom des gènes", x = "Nombre de variants", y = "VCF") +
theme(panel.grid = element_blank(),plot.title = element_text(hjust = 0.5)) +
labs(title = "Nombre de variants par gène et par VCF")
ggplotly(plot2)Le nombre de variants est identique dans chaque VCF pour les gènes UTS2 et LPL mais le nombre varie pour le plus grand gène UPHA5.
Nombre et type de variants par gène et par VCF
nombre_variants_par_gene$vcf_nom_genes = paste(nombre_variants_par_gene$vcf, nombre_variants_par_gene$gene, sep="_")
plot3 = ggplot(nombre_variants_par_gene, aes(fill=type, y=vcf_nom_genes, x=nombre)) +
geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
scale_fill_manual(values = c("#bebada", "#ffffb3", "#8dd3c7", "#170f0f")) +
labs(title = "Nombre de variants par gène, par VCF et par type de variant") +
labs(fill = "Type de variants", x = "Nombre de variants", y = "VCF") +
theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
ggplotly(plot3)La répartition des types de variants par gène et par VCF est concordante. Aucune insertion dans UTS2 dans les 4 VCF.
Exploration des SNV
Nombre et type de SNV par VCF
plot4 = ggplot(nombre_type_snv[!is.na(nombre_type_snv$type), ], aes(fill = type, y = vcf, x = nombre))+
geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
scale_fill_manual(values = c("#fc9272", "#9ebcda")) +
labs(title = "Type de SNV par VCF", x = "Type de substitutions", y = "Nombre de substitutions") +
labs(fill = "Type de SNV") +
theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
plot4
En génétique humaine, il y a un environ 2 fois plus de transitions que de transversions. Sur le plot, on voit environ 2/3 de transitions et 1/3 de transversions ce qui est cohérent. Valeurs cohérentes dans tous les VCF.
Nombre et type de transitions/transversions par VCF et par gène
plot5 = ggplot(nombre_type_snv[!is.na(nombre_type_snv$type_snv), ], aes(fill = type_snv, y = vcf_gene, x = nombre)) +
geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
labs(title = "Répartition des types de SNV par VCF et par gène", x = "Nombre", y = "VCF et gène") +
scale_fill_manual(values = c("#a6cee3", "#1f78b4","#b2df8a", "#33a02c", "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a", "#ffff99", "#b15928")) +
labs(fill = "Type de SNV", x = "Nombre", y = "VCF") +
theme(plot.title = element_text(hjust = 0.5, size=14, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
theme(panel.grid = element_blank())
ggplotly(plot5)Une répartition des types de SNV cohérente entre VCF pour chaque gène.
Exploration des délétions/insertions
plot6 = ggplot(vcf_all %>% filter (type == "insertion" | type == "deletion"), aes(x = taille_delins, y = vcf)) +
geom_boxplot(aes(fill=type)) +
theme_classic() +
labs(x = "Taille des delins", y = "VCF", title = "Boxplots de la taille des insertions et délétions par VCF") +
theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
labs(fill = "Type de delins") +
scale_fill_manual(values = c("deletion" = "#bebada", "insertion" ="#ffffb3"))
plot6
La distribution des tailles des délétions est équivalente entre les 4 VCF. Les VCF 12x, 24x et 40x contiennent plus de valeurs aberrantes pour les insertions et délétions que le VCF de référence. La taille maximale des insertions du VCF de référence est plus petite que pour les 3 VCF obtenus avec le pipeline. L’interquartile est plus petit également pour la taille des insertions du VCF de référence ce qui signifie que les valeurs des insertions sont plus regroupées autour de la médiane. Il y a moins de disparité des valeurs pour le VCF de référence.
plot7 = ggplot(vcf_all %>% filter (type == "insertion"), aes(x = taille_delins, y = vcf)) +
geom_boxplot(aes(fill=gene)) +
theme_classic() +
labs(x = "Taille des insertions", y = "VCF", title = "Boxplots de la taille des insertions par VCF et par gène") +
theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
labs(fill = "Noms des gènes") +
scale_fill_manual(values = c("EPHA5" = "#fdc086", "LPL" ="#91bfdb", "UTS2" = "#b2df8a"))
plot8 = ggplot(vcf_all %>% filter (type == "deletion"), aes(x = taille_delins, y = vcf)) +
geom_boxplot(aes(fill=gene)) +
theme_classic() +
labs(x = "Taille des délétions", y = "VCF", title = "Boxplots de la taille des délétions par VCF et par gène") +
theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
labs(fill = "Noms des gènes") +
scale_fill_manual(values = c("EPHA5" = "#fdc086", "LPL" ="#91bfdb", "UTS2" = "#b2df8a"))
plot_grid(plot7, plot8)
Exploration de la profondeur de lecture des variants
plot9 = ggplot(vcf_12x_24x_40x, aes(x = as.numeric(vcf_12x_24x_40x$DP), y = vcf)) + geom_boxplot(aes(fill=gene)) + theme_classic() +
theme_classic() +
labs(x = "Profondeur de lecture des variants", y = "VCF", title = "Boxplots de la profondeur de lecture des variants par gène") +
theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
labs(fill = "Noms des gènes") +
scale_fill_manual(values = c("UTS2" = "#b2df8a", "LPL" ="#91bfdb", "EPHA5" = "#fdc086"))
plot9## Warning: Use of `vcf_12x_24x_40x$DP` is discouraged.
## ℹ Use `DP` instead.
plot10 = ggplot(dp_12x, aes(x = DP, y = nombre)) +
geom_histogram(stat = "identity", fill = "#a6cee3") +
geom_vline(xintercept = mean(dp_12x$DP), color = "#fd8d3c", linetype = "dashed", size = 1) +
labs(x = "Profondeur de lecture", y = "Nombre de variants", title = "VCF 12x") +
theme(axis.text.x = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5, size = 14, face = "bold")) +
theme(panel.grid = element_blank())## Warning in geom_histogram(stat = "identity", fill = "#a6cee3"): Ignoring
## unknown parameters: `binwidth`, `bins`, and `pad`## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.plot11 = ggplot(dp_24x, aes(x = DP, y = nombre)) +
geom_histogram(stat = "identity", fill = "#1d91c0") +
geom_vline(xintercept = mean(dp_24x$DP), color = "#fd8d3c", linetype = "dashed", size = 1) +
labs(x = "Profondeur de lecture", y = "Nombre de variants", title = "VCF 24x") +
theme(axis.text.x = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5, size = 14, face = "bold")) +
theme(panel.grid = element_blank())## Warning in geom_histogram(stat = "identity", fill = "#1d91c0"): Ignoring
## unknown parameters: `binwidth`, `bins`, and `pad`plot12 = ggplot(dp_40x, aes(x = DP, y = nombre)) +
geom_histogram(stat = "identity", fill = "#253494") +
geom_vline(xintercept = mean(dp_40x$DP), color = "#fd8d3c", linetype = "dashed", size = 1) +
labs(x = "Profondeur de lecture", y = "Nombre de variants", title = "VCF 40x") +
theme(axis.text.x = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5) ) +
theme(panel.grid = element_blank())## Warning in geom_histogram(stat = "identity", fill = "#253494"): Ignoring
## unknown parameters: `binwidth`, `bins`, and `pad`plot_grid(plot10, plot11, plot12, ncol = 3, nrow = 1)
Exploration de la VAF des VCF 12x, 24x et 40x
On ne fera pas d’exploration de la VAF ni de la profondeur DP pour le VCF de référence car il n’y a pas d’explication aux valeurs de DP, AD, ADALL dans le VCF. Les valeurs ne sont pas cohérentes avec ce qui est observé dans le bam de référence associé au VCF dans IGV.
VAF minimale détectée dans chaque VCF
cat("VAF minimale dans VCF 12x :", min(vaf_12x$vaf), "%", "\n")## VAF minimale dans VCF 12x : 16.7 %cat("VAF minimale dans VCF 24x :", min(vaf_24x$vaf), "%", "\n")## VAF minimale dans VCF 24x : 14.3 %cat("VAF minimale dans VCF 40x :", min(vaf_40x$vaf), "%", "\n")## VAF minimale dans VCF 40x : 16.3 %plot13 = ggplot(vcf_vaf, aes(x = nombre_variants, y = cat_vaf , fill = vcf)) +
geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
labs(title = "Répartition du nombre de variants par catégorie de VAF", x = "Nombre de variants", y = "Catégorie de VAF") +
scale_fill_manual(values = c("#a6cee3", "#1d91c0", "#253494")) +
labs(fill = "VCF") +
theme(plot.title = element_text(hjust = 0.5, size=14, face="bold"), plot.subtitle = element_text(hjust = 0.5)) +
theme(panel.grid = element_blank()) +
theme_classic()
ggplotly(plot13)Exploration des variants en commun dans les 4 VCF
variants_list = list(vcf_12x = unique(vcf_12x$variants), vcf_24x = unique(vcf_24x$variants), vcf_40x = unique(vcf_40x$variants), vcf_ref = unique(vcf_ref$variants))
plot14 = upset(fromList(variants_list), order.by = "freq", text.scale = 2, point.size = 5, sets.bar.color = "#999999", main.bar.color = "#fc8d62")
plot14
plot15 = display_venn(variants_list, category.names = c("vcf_12x" , "vcf_24x" , "vcf_40x", "vcf_ref"),
fill = c("#a6cee3", "#1d91c0", "#253494", "#E69F00"), cat.cex = 1.2,
cat.fontface = "bold",
cat.default.pos = "outer",
cat.dist = c(0.1, 0.1, 0.1, 0.1))
plot16 = ggplot(cat_vcf_all_nombre, aes(fill=type, x=nombre, y=categorie)) +
geom_bar(position = "stack", stat="identity") +
scale_fill_manual(values = c("#bebada", "#ffffb3", "#8dd3c7", "#170f0f")) +
labs(title = "Nombre de variants par catégories et par type de variants") +
labs(fill = "Type de variants", x = "Nombre de variants", y = "Catégories de VCF") +
theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))plot17 = ggplot(cat_vcf_all_genes, aes(fill=gene, x=nombre, y=categorie)) +
geom_bar(position = "stack", stat="identity") +
scale_fill_manual(values = c("UTS2" = "#b2df8a", "LPL" ="#91bfdb", "EPHA5" = "#fdc086")) +
labs(title = "Nombre de variants par catégories et par gènes") +
labs(fill = "Type de variants", x = "Nombre de variants", y = "Noms des gènes") +
theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
plot_grid(plot16, plot17, ncol = 1, nrow = 2)
Exploration des variants en commun dans les 3 VCF générés par le pipeline
variants_list_12x_24x_40x = list(vcf_12x = unique(vcf_12x$variants), vcf_24x = unique(vcf_24x$variants), vcf_40x = unique(vcf_40x$variants))
plot18 = display_venn(variants_list_12x_24x_40x, category.names = c("vcf_12x" , "vcf_40x", "vcf_24x"),
fill = c("#a6cee3", "#1d91c0", "#253494"), cat.cex = 1.2,
cat.fontface = "bold",
cat.default.pos = "outer",
cat.dist = c(0.1, 0.1, 0.1))
Exploration des métriques de qualité: datas générés par Hap.py d’Illumina
recall (sensibilité), précision (spécificité) et score F1. Ce sont des métriques qui permettent de voir si un modèle est performant. Le score F1 combien recall et précision Plus le recall est élevé, plus le modèle repère des positifs. Plus la précision est élevée, moins le modèle se trompe sur les positifs (moins il y a de faux positifs). Plus le score F1 est élevé, plus le modèle est performant.
Courbe ROC
happy_12x_summury = read.table("/home/elodie/Documents/Analysis/Happy_12x/Happy_12x/ref-12x.summary.csv", header=TRUE, sep=',')
happy_24x_summury = read.table("/home/elodie/Documents/Analysis/Happy_24x/Happy_24x/ref-24x.summary.csv", header=TRUE, sep=',')
happy_40x_summury = read.table("/home/elodie/Documents/Analysis/Happy_40x/Happy_40x/ref-40x.summary.csv", header=TRUE, sep=',')plot_roc_snp = plot_data(roc_snp, is.PR=TRUE)## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.plot_roc_indel = plot_data(roc_indel, is.PR=TRUE)
plot_grid(plot_roc_snp, plot_roc_indel, nrow=2, ncol=1)